A Direct Elliptic Solver Based on Hierarchically Low-rank Schur Complements

Cyclic reduction was conceived for the solution of tridiagonal linear systems, such as the one-dimensional Poisson equation [11]. Generalized to higher dimensions, it is known as block cyclic reduction (BCR) [4]. It can be used for general (block) Toeplitz and (block) tridiagonal linear systems; however, it is not competitive for large problems, because its arithmetic complexity grows superlinearly. Cyclic reduction can be thought of as a direct Gaussian elimination that recursively computes the Schur complement of half of the system. Schur complement computations have complexity that grows with the cost of the inverse, but by considering a tridiagonal system and an even/odd ordering, cyclic reduction can decouple the system in such a way that the inverse of a large block is the block-wise inverse of a collection of independent smaller blocks. This addresses the most expensive step of the Schur complement computation in terms of operation complexity and does so in a way that launches concurrent subproblems. Its concurrency feature in the form of recursive bisection makes it interesting for parallel environments, provided that its arithmetic complexity can be improved. We address the time and memory complexity growth of the traditional cyclic reduction algorithm by approximating dense blocks as they arise with hierarchical matrices (H-matrices). The effectiveness of the block approximation relies on the rank structure of the original matrix. Many relevant operators are known to have low rank off the diagonal. This philosophy follows recent work discussed below, but to our knowledge this is the first demonstration of the utility of complexity-reducing hierarchical substitution in the context of cyclic reduction.